/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

// forceNbdES.c -- functions to evaluate (nonbonded) electrostatic pair interactions

#include "definitions.h"

// Electric conversion factor relating mechanical to electrical quantities -
// see also Gromacs user's manual, "Definitions and Units" chapter
static const real ONE_OVER_4_PI_EPSILON0_IN_BRAHMS_UNITS =  1 / ( 4 * PI * EPSILON0_IN_BRAHMS_UNITS ); 

// Coulomb interaction, Shifed-Force truncation
void CalcChgChgInteraction_SF( real *ijEnergy, // out - potential energy
			       VecR *ijForce,  // out - force vector
			       real prodCharges, // in - product of charges (iCharge*jCharge)
			       const real rCut, // in - cutoff distance
			       const real rr,  // in - intersite distance squared
			       const VecR rij ) // in - intersite distance vector 
{
  real minFacSq, minusFac, plusFac, r, uElecij;
  
  prodCharges *= ONE_OVER_4_PI_EPSILON0_IN_BRAHMS_UNITS; // scaling according to the electric conversion factor
  r = sqrt( rr ); 
  minusFac = 1. - r / rCut;
  plusFac =  1. + r / rCut;
  minFacSq = Sqr( minusFac );
  uElecij = prodCharges / r;
  *ijEnergy = uElecij * minFacSq;
  VSCopy( *ijForce, uElecij * minusFac * plusFac / rr, rij ); 
}

// Charge-dipole interaction, Shifed-Force truncation
void CalcChgDipInteraction_SF( real *ijEnergy, // out - potential energy
			       VecR *forceChgDip,// out - force on charge
			       VecR *torqueDip,   // out - torque on dipole
			       real prodMagChgDip, // in - product of magnitudes
			       const VecR e, /* in - dipole orientation vector */
			       const real rCut, // in - cutoff distance
			       const VecR rChgDip, // in - distance (rChg-rDip)  
			       const real rr )    // in - rChgDip squared
			     
{
  real elMagOvrrr, er, rrr, rrrC, rOvrC_Cube, rOvrC_Sq, sfFac;

  prodMagChgDip *= ONE_OVER_4_PI_EPSILON0_IN_BRAHMS_UNITS; // scaling according to the electric conversion factor
  rrr = sqrt( rr ) * rr;
  rrrC = rCut * Sqr( rCut ); 
  rOvrC_Sq = rr / Sqr( rCut );
  rOvrC_Cube = rrr / rrrC;
  elMagOvrrr = prodMagChgDip / rrr;
  sfFac = 1. - 3. * rOvrC_Sq + 2. * rOvrC_Cube;
  er = VDot( e, rChgDip );
  *ijEnergy = elMagOvrrr * sfFac * er;
  VSSAdd( *forceChgDip, 3. * er / rr * ( 1. - rOvrC_Sq ), rChgDip, -sfFac, e ); 
  VScale( *forceChgDip, elMagOvrrr ); 
  VCross( *torqueDip, rChgDip, e ); 
  VScale( *torqueDip, elMagOvrrr * sfFac );
}

/* Dipolar interaction [allen & tildesley, appendix C] with cubic switching for rSwitch < r < rCut */
void CalcDipDipInteraction_CS( real *ijEnergy, // out - potential energy
			       VecR *ijForce,  // out - force 
			       VecR *ijTorque, // out - torque
			       VecR *jiTorque, // out - torque
			       real dipMagProd, // in - product of dipole magnitudes
			       const VecR iOrient, // in - "i" orientation (unit) vector 
			       const VecR jOrient, // in - "j" orientation (unit) vector 
			       const real rCut, // in - cutoff distance
			       const VecR rij, // in - intersite distance vector 
			       const real rr ) // in - intersite distance squared
{
  real cosGamma, cosTheta_i, cosTheta_j, CrUmrL, r, rrr, rUmr, s, sFactor;
  VecR cosVec1, cosVec2, cosVec, dipTrq1, dipTrq2, rCosVec, rijUnit, switchGradient;

  const real rSwitch = 0.9 * rCut; // [Leach, Molecular modelling, p.332, 2nd ed]
  
  dipMagProd *= ONE_OVER_4_PI_EPSILON0_IN_BRAHMS_UNITS; // scaling according to the electric conversion factor
  r = sqrt ( rr ); 
  rrr = rr * r;
  VSCopy( rijUnit, 1./ r, rij );
  cosGamma = VDot( iOrient, jOrient ); /* C.15 */
  cosTheta_i = VDot( iOrient, rijUnit );
  cosTheta_j = VDot( jOrient, rijUnit );
  *ijEnergy = dipMagProd * ( cosGamma - 3. * cosTheta_i * cosTheta_j ) / rrr; /* (C.16) */
  /* dipolar force calculation  - from C.26 a, allen & tildesley */
  VSCopy( cosVec1, cosTheta_j, iOrient );
  VSCopy( cosVec2, cosTheta_i, jOrient );
  VAdd( cosVec, cosVec1, cosVec2 );
  VSCopy( rCosVec, ( cosGamma - 5. * cosTheta_i * cosTheta_j ), rijUnit );
  VAdd( *ijForce, rCosVec, cosVec );
  VScale( *ijForce, 3. * dipMagProd / Sqr( rr ) );
  /* dipolar torque calculation  - from C.26 b-c, allen & tildesley */
  VCross( dipTrq1, iOrient, jOrient );
  VCross( dipTrq2, iOrient, rij );
  VScale( dipTrq2, 3. *  cosTheta_j / r );
  VSub( *ijTorque, dipTrq2, dipTrq1 );
  VScale( *ijTorque, dipMagProd / rrr );
  VCross( dipTrq1, jOrient, iOrient );
  VCross( dipTrq2, jOrient, rij ); /* doubt: why not rji? */
  VScale( dipTrq2, 3. * cosTheta_i / r );
  VSub( *jiTorque, dipTrq2, dipTrq1 );
  VScale( *jiTorque, dipMagProd / rrr );
  if ( r > rSwitch ) { /* switching for rSwitch < r < rCut*/
    rUmr = rCut - r;
    CrUmrL = Cube( rCut - rSwitch );
    sFactor =  6. * rUmr * ( rSwitch - r ) / ( CrUmrL * r );
    VSCopy( switchGradient, sFactor, rij );
    s = Sqr( rUmr ) * ( rCut + 2. * r - 3. * rSwitch ) / CrUmrL;
    VScale( *ijForce, s );
    VScale( *ijTorque, s );
    VScale( *jiTorque, s );
    VVSAdd( *ijForce, -1. * *ijEnergy, switchGradient );
    *ijEnergy *= s;
  }
}
